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1.  Introduction 


Several  signal  detection  and  estimation  problems  arise  in  monitoring  teleseismic  and  re¬ 
gional  events  using  data  from  long  period  and  short  period  arrays.  Both  on-Mne  and  off-line 
detection  and  estimation  software  must  be  able  to  operate  successfully  in  the  presence  of 
either  relatively  high  noise  or  interfering  events.  Besides  lowering  detection  thresholds,  the 
interfering  events  will  produce  distortions  in  velocity  and  azimuth  estimators  resulting  in 
poor  locations  for  teleseismic  events  and  in  an  inability  to  sort  out  velocities  for  the  differ¬ 
ent  phases  for  near  and  far  regional  arrivals  on  the  short  period  arrays.  It  should  also  be 
noted  that  amplitude  or  power  measurements  made  for  discrimination  purposes  will  also  be 
distorted. 

Fitting  data  from  an  array  of  distributed  sensors  to  a  multiple-signal  model  for  the  purpose  of 
extracting  possible  mixed  signals  is  a  method  used  routinely  in  seismic  data  processing,  e.g., 
Smart(1972).  In  general,  however,  such  procedures  do  not  constitute  detectors.  They  have 
served,  among  other  purposes,  to  disentangle  mixed  (simultaneously  arriving)  long-period 
surface- wave  signals  when  one  of  them  was  sought  as  a  critical  discriminant  to  distinguish  an 
explosion  from  a  possible  earthquake.  Such  processors  have  had  no  independent  detection 
statistics,  and,  consequently,  for  verification  and  validation  of  waveforms  extracted  by  their 
use,  alternative  means  are  employed.  Those  means  include  the  search,  in  the  same  time 
interval,  for  body- waves  of  passible  associated  seismic  events  in  the  records  of  a  complemen¬ 
tary  network  of  short-period  seismic  stations.  Any  tentatively  identified  surface  waveform, 
recovered  by  multiple-signal  modehng,  may  be  verified  by  association  with  one  of  the  seismic 
events  detected  and  located  by  that  search.  If  a  surface  waveform  so  extracted  meets  the 
criteria  of  arrival  time,  dispersion,  and  magnitude  predicted  from  the  estimated  location, 
magnitude  and  time  of  occurrence  of  a  seismic  event  discovered  in  the  short-period  data,  the 
surface-wave  is  associated  with  that  event  and  treated  as  'genuine'- 

To  illustrate,  consider  Figure  1,  which  contains  a  verifiable  mixture  of  signals  from  two 
earthquakes,  one  from  the  south  of  Africa  and  the  other  from  the  Philippines,  observed  at 
a  long  period  array  in  Korea.  As  we  will  note  later,  conventional  methods  applied  to  this 
data  indicate  a  single  event  from  somewhere  between  the  two  observed  events.  Rather  than 
discovering  such  incorrect  results  using  an  alternate  network  or  short-period  records,  it  would 
be  strongly  preferable  to  have  a  detection  statistic  that  sorted  out  multiple  arrivals  from  the 
long-period  records  in  Figure  1.  -We  show  that  the  sequential  detection  procedure  given  in 
this  paper  identifies  the  two  signals,  as  well  as  a  complementary  unidentified  component  that 
may  be  another  event  or  noise. 

Another  powerful  advantage  of  the  multiple-signal  detector  is  its  immunity  to  interference 
from  persistent  coherent  noise,  natural  or  anthropogenic,  from  such  phenomena  as,  say,  the 
microseisms  of  storms  at  sea,  or  from  power  plant  emanations.  Conventionally  persistent, 
nuisance  coherences  are  dealt  with  by  notch  filters  to  prevent  the  detector  from  responding 
to  the  coherent  noise,  and  from  missing  smaller  signals  in  adjoining  frequencies.  These  notch 
filters  generally  operate  right  in  the  midst  of  the  band  of  frequencies  of  interest,  and  that 
critical  portion  of  the  band  is  sacrificed  to  accommodate  the  vulnerabilities  of  the  detector 
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Figure  1:  Mixture  of  signals  from  two  earthquakes  from  south  of  Africa 
and  the  Philippines  observed  at  USAEDS  long-period  seismic 
array  in  Korea.  Correct  geodesic  back-azimuths  are  226  degrees 
and  198  degrees. 


(Clauter,  2004).  Finally,  we  note  another  benefit  of  the  sequential  F-detector  proposed  in 
this  paper  is  its  use  of  the  noise  computed  during  the  signal  window.  Conventional  detection 
using  the  ratio  of  the  short-term  average  mean-square  error  (STA)  to  the  corresponding 
long-term  average  (LTA)  requires  that  the  LTA  be  refreshed  after  the  onset  of  a  signal  to 
remain  relevant,  causing  a  blind  window  where  the  detector  will  miss  valid  signals. 
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Approaches  to  detecting  signals  on  arrays  all  focus  on  the  basic  model  that  expresses  the 
observed  channel  as  sums  of  delayed  signals  and  a  unique  noise  process.  The  delays  are 
functionally  dependent  on  velocity  and  azimuth  if  the  signals  are  propagating  plane  waves, 
and  this  is  the  assumption  that  is  usually  made.  Methods  that  are  commonly  in  use  for  ana¬ 
lyzing  such  data,  when  a  single  signal  is  assumed  to  be  present,  can  be  roughly  categorized  as 
(a)  beam-forming  and  plotting  the  power  as  a  function  of  slowness,  which  can  be  converted 
to  estimators  of  velocity  and  azimuth, (b)  Capon’s  estimator  (see  Capon,  1969,  Capon  and 
Goodman,  1970),  (c)  beam-forming  converted  to  an  F-statistic  by  dividing  by  an  estimator 
of  the  noise  power  (see  Melton  and  Bailey,  1957,  Blandford,  1970,  Shumway,  1970,  1971, 
1983,  2000),  (d)  Multiple  Signal  Characteristic  (MUSIC) (Schmidt,  1979,  Stoica  and  Neho- 
rai,  1989),  and  (e)  cross  correlation  (TVibuleac  and  Herrin,  1997).  Several  algorithms  such 
as  the  sequential  F-detector  proposed  here  and  the  multiple  signal  characteristic  (MUSIC) 
algorithm  are  available  that  oflFer  promise  for  handling  array  data  with  low  signal-to-noise 
ratios  and  contamination  from  interfering  signals.  In  this  proposal,  we  exhibit  the  perfor¬ 
mance  of  currently  available  algorithms  on  teleseismic  and  regional  data  containing  mixed 
signals  and  demonstrate  the  superior  performance  of  the  sequential  F-statistic.  A  sequen¬ 
tial  analysis  of  power  using  the  F-statistic  is  proposed  that  estimates  the  correct  number  of 
signals  and  their  velocities  and  azimuths.  This  is  contrasted  with  results  using  conventional 
f-k  estimators  that  do  not  handle  the  mixed  signal  case. 

2.  The  Multiple  Signal  Model 

The  usual  model  expresses  the  received  data  =  1, ...  ,N,t  =  I ...  ,n  a,t  N  sensors 

M 

yj{i)  =  Jlsk{t  +  r'fik)  +Vj{t)  ( 1 ) 

fc=i 

as  the  sum  of  M  propagating  plane  waves  with  time  delays  7)  (Ok)  =  for  the  signal 
at  the  sensor  where  rj  is  the  two  dimensional  coordinate  of  sensor  j  in  km  and  0k  is  the 
two-dimensional  slowness  vector  of  signal  k  in  sec/km.  The  two-dimensional  slowness  vector 
—  (^ifci  ^2fc)^  is  related  to  the  velocity  14  =  ||^fei|“^  and  azimuth  a*  =  tan~^(02(t/^ife))- 

Because  the  plane  wave  model  formulates  more  easOy  in  the  frequency  domain,  we  may 
consider  a  model  for  the  discrete  Fourier  transforms,  say 

M 

exp{27rwuf  rj0k}Sk{ujij  +  (2) 

fc=i 

evaluated  at  ^  =  1.  ..L  frequencies  of  the  form  ==  l/n  m  the  neighborhood  of  some 
target  frequency  tu,  measured  in  cycle  per  point.  A  convenient  representation  of  the  model 
is  obtained  by  stacking  the  sensors  in  an  AT  x  1  vector  Yi  =  {Yi{uJe)  ...,Y  Then,  the 

above  equation  becomes 

Ye  =  Zeie)Se  +  Ve  (3) 

where  Se  —  (Siicoi), . . . ,  denotes  the  unknown  vector  of  signals  at  frequency  and 

Ze(0)  =  |exp  {27riaj{  rj0k},j  =  1,. . .  ,N,k  ~  1 . . .  ,Mj  (4) 
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is  an  iV  X  M  matrix  that  defines  the  way  the  plane  waves  map  into  the  observed  sensor 
elements.  Suggested  detectors  are  based  on  the  statistical  structure  assumed  for  the  elements 
of  the  model  (2)  under  the  fixed  and  stochastic  signal  assumptions. 

Suppose  that  the  signal  is  deterministic  and  the  noise  has  N  x  N  spectral  matrix  /i,(w). 
Then,  it  is  clear  that  asymptotically,  the  inferences  can  be  based  on  a  model  that  assumes 
that  Ye,  i  =  1, . . . ,  L  axe  independently  distributed  with  mean  Z{Q)Se  and  covariance  matrix 
fv  Various  specializing  assumptions  lead  to  the  classic  beam  and  to  the  various  forms  of  the 
F-detector.  We  remark  that  it  is  clear  that  the  deterministic  model  has  the  nonlinearity  in 
the  mean  function.  The  primary  interest  will  be  in  the  parameter  matrix  0  =  (0i ,  02,  ■  ■  ■ ,  0m) 
which  enters  non-linearly  in  the  mean  function, 

EYe)  =  Ze(0)Se. 

If  the  signal  is  random  with  qxq  spectral  matrix  /s(a?),  Ye,C  =  l,. . .  ,L  will  be  independently 
distributed  with  mean  zero  and  covariance  matrix 

E,(a;)  Ze{e)Uuj)Zm  +  (5) 

and  the  nonlinearity  in  the  slowness  parameters  is  concentrated  in  the  covariance  structure. 
Such  a  model  can  be  used  to  argue  in  the  single-signal  case  for  Capon’s  (1969)  high-resolution 
estimator  and  for  the  MUSIC  estimator  of  Schmidt  (1979).  Furthermore,  a  maximum  like¬ 
lihood  estimator  can  be  derived  by  maximizing  the  quasi-Gaussian  log  likelihood  implied 
under  the  complex  Gaussian  assxunption  (see  Shumway  et  al,  1999). 

There  are  a  number  of  conventional  approaches  to  estimating  the  slowness  parameters  and 
hence,  the  derived  velocities  and  azimuths.  Most  common  are  the  beam  power  and  F- 
statistics  described  in  Section  3.1  below  and  variations  on  the  correlation  method,  described 
in  (3.4).  Note  that  the  F-statistic  is  a  monotone  function  of  the  statistic,  semblence.  We  also 
adapt  the  Multiple  Signal  Characteristic  (MUSIC)  algorithm  of  Schmidt  (1979)  to  seismic 
arrays.  Finally,  we  look  at  the  Capon  (1969)  "high-resolution"  estimator  as  another  variation 
based  on  the  eigen  vectors  and  eigen  values  of  the  spectral  matrix. 

The  following  Section  4  introduces  the  multiple  signal  approach  using  the  partial  F-statistic 
which  is  shown  to  improve  on  the  conventional  astimators  in  Section  5. 

3.  Conventional  Approaches 

We  summarize  first  some  common  approaches  to  determining  the  velocity  and  azimuth  of 
a  possible  mixture  of  propagating  signals.  The  model  is  expressed  in  terms  of  slowness  in 
seconds  per  km  0i  —  (0l,02)^  where  6i  and  02  are  the  horizontal  and  vertical  components. 
Velocity  c  =  1/||0||  and  back  azimuth  a  —  tan''^(02/0i)  are  the  propogating  plane  wave 
velocities  in  km/sec  and  back  azimuths  from  due  north.  We  will  see  later  that  the  techniques 
described  below  do  not  work  well  on  some  common  mixtures. 
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3.1  Beam  Power  and  F-statistics 

Consider  the  single  signal  version  of  (1),  namely,  take  Af  =  1  so  that  the  model  becomes 


Vj  (t)  =  >51  (t  +  r'jOi )  +  Vj  (t)  (6) 

for  the  sensor.  Transforming  this  model  to  the  frequency  domain  and  stacking  the 
observed  series  yields 

Yi  =  zn(ei)Sn+Ve  (7) 

where 

ZiiiOi)  ^  {exp{2mu}e  YiOi),  •  •  • , exp{27ria;^  (8) 

is  the  frequency  domain  version  of  the  beam-steering  vector  when  it  is  concentrated  at  the 
true  slowness  6i.  The  beam  probe 


(9) 


will  be  large  eLt0  —  0i  and  will  be  small  off  the  beam.  If  the  spectral  matrix  /„(w)  =  Sy{u)Iff 
is  diagonal  with  equal  spectra  on  all  channels,  then  the  beam  power 


Pam  =  E  lamp 


f=i 


(10) 


is  distributed  proportionally  to  a  chi-squared  random  variable  with  2L  degrees  of  freedom 
when  0  =  01  the  true  slowness  value.  Unfortunately  the  proportionality  constant  is  a  function 
of  the  unknown  noise  spectrum  The  noise  spectrum  can  be  estimated  from  a  sample 

of  data  prior  to  the  signal,  but  the  problems  with  that  practice  have  already  been  mentioned. 
Note  that  L  =  DT  where  B  is  the  bandwidth  in  Hz,  and  T  is  the  sample  length  in  seconds. 

It  is  also  the  case  that  the  value  of  0  that  maximizes  the  beam-power  (6)  also  minimizes  the 
mean  squared  error 

SSE{0)  =  E \\Ye  -  zn{0)De{0)\\\  (11) 

t=i 


and  it  can  be  shown  that  the  hkelihood  ratio  test  of  no  signal  in  the  model  (5)  leads  to  the 
F-Statistic 


F{0)  =  c 


SSE{0)  ’ 


(12) 


which  converges  to  anon-central  F-distribution  with  dfi  —  2(T-fl)  and  d/2  =  2[L(iV  — 1)  — 1] 
degrees  of  freedom  as  0  »  0i.  The  constant  is  c  =  d/2/dfi  The  advantage  of  the  F-statistic 
is  that  its  distribution  does  not  depend  on  nuisance  parameters  so  that  thresholds  can  be  set 
directly  from  the  tabulated  values.  Furthermore,  an  estimator  of  the  noise  power  is  available 
from  the  signal  window  as  P„(a;)  =  SSE{0)/N. 
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3.2  Capon’s  High  Resolution  Method 

Estimators  available  in  the  engineering  literature  depend  primarily  on  the  sample  spec¬ 
tral  matrix  and  on  the  average  wave-number  value  over  a  frequency  band.  For  example 
Capon’s(1971)  estimator  uses  the  sample  spectral  matrix 

s,=Env;  (13) 

as  input  to  the  proposed  estimator  obtained  by  maximizing  the  function 

C{e)  =  (14) 

which  has  a  distribution  proportional  to  the  chi-squared  distribution  with  2{L  —  +  1) 

degrees  of  freedom  (see  Capon  and  Goodman,  1970).  The  test  statistic  is  evaluated  at  at  an 
integer  near  the  mean  frequency,  say  1.  The  distribution  here  also  depends  on  the  unknown 
signal  and  noise  spectra  and  P„,  or  on  the  noise  spectrum  only  when  there  is  no  signal, 
i.e.,  Pg  =  0 

3.3  The  MUSIC  Algorithm 


The  MUSIC  approach  suggested  by  Schmidt  (1979)  has  been  extended  to  arrays  by  Shumway 
(2002)  for  possible  use  in  detecting  infrasound  signals.  The  approach  is  based  on  orthogonal¬ 
ity  properties  of  the  eigen  vectors  of  the  spectral  matrix  under  a  simple  multivariate  random 
signal  model.  In  general,  the  test  statistic  (11)  is  expected  to  have  large  values  under  a 
model  that  assumes  an  appropriate  number  of  input  signals. 


An  estimator  based  on  the  eigen  vectors  of  Sy  forms  the  main  component  of  the  Multiple 
Signal  Characteristic  (MUSIC)  detector  proposed  by  Schmidt  (1979).  A  good  summary  of 
the  statistical  properties  of  this  estimator  is  in  Stoica  and  Nehorai  (1989).  In  particular, 
letting  ei,e2, . . .  jCat  be  the  eigen  vectors  of  Ey,  take  the  MUSIC  detector  as  the  maximizer 


of 


M{e)  =  \z\{9){  f:  eie*\zi{9) 


•j=M+l 


-1 


(15) 


where  we  assume  that  there  are  M  possible  signals. 


The  approach  to  determining  the  number  of  signals  and  their  velocities  and  azimuths  reduces 
to  trying  various  values  of  M  in  the  above  equation,  say  M  =  1, 2, 3, 4,  looking  for  peaks  in 
the  plotted  function  M{9). 

It  is  interesting  to  note  that  the  Capon  detectors  can  also  be  written  in  terms  of  the  eigen 
values  and  eigen  vectors  of  the  spectral  matrix.  We  note  that  (14)  becomes 


C{9)  = 


The  equation  exhibits  the  estimator  in  terms  of  a  matching  of  the  probe  vector  with  the 
eigen  vectors. 
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3.4  Correlation  Methods 


There  are  several  simple  variations  that  use  the  cross  correlation  pairs  measured  over  the 
array  to  estimate  the  velocity  and  azimuth.  Note  first  that  the  single  signal  model  (5)  implies 
that  the  sample  cross  correlation  function  between  yj{t)  and  yk{t)  will  tend  to  be  maximized 
at  lag  (r^  -  rfeX^i. 

One  simple  method  might  be  to  simply  average  the  cross  correlation  functions  which  turns 
out  to  be  nearly  the  same  as  broad  band  beam-forming.  A  somewhat  better  method  is  to 
think  of  the  estimated  lag  Tjk  that  maximizes  the  cross  correlation  between  yj{t)  and  yk{t) 
to  be  represented  by  the  model 


Tjk  =  (rj  -  rkYOi  +  Cjk  (16) 

for  j  ^  k  —  1 ...  ,N  and  minimize  some  objective  function  of  the  residuals  ejk.  Tribuleac 
and  Herrin  (1997)  use  the  sum  of  the  absolute  errors,  robustified  by  dropping  out  some  the 
largest  ones. 

The  computations  done  later  here  simply  minimized  the  sum  of  squared  errors  with  no 
robustification.  For  example,  just  chose  6>i  to  minimize 

sse  YliTjh  -  (rj  -  rkYOif.  (17) 

j<k 

One  problem  with  these  approaches  is  that  standard  errors  of  the  estimator  for  0i  are  not 
easy  to  come  by.  The  distribution  of  the  maximizing  values  Tjk  is  not  available,  particularly 
when  it  is  noted  that  the  errors  should  be  highly  correlated  because  of  common  sensors  used 
in  calculating  the  cross  correlations,  for  example,  puiPiz,  •  •  • ,  Pin  all  involve  the  first  sensor. 


3.5  Remarks 

Given  the  erratic  performance  of  conventional  estimators  to  be  demonstrated  later  in  the 
case  of  mixed  signals,  developing  a  means  for  automatically  detecting  and  isolating  multiple 
interfering  signals  should  have  high  priority.  First  of  all,  one  would  like  to  have  a  reliable 
means  of  detecting  a  single  signal  over  noise,  and  all  methods  appear  to  do  an  acceptable  job 
of  indicating  local  maximums  that  agree.  The  question  of  statistical  significance  of  the  single¬ 
signal  detection  rests  on  the  probability  of  obtaining  a  larger  value  for  any  one  of  detectors 
(a)-(d)  or  the  cross  correlation  method.  A  satisfactory  solution  to  the  problem  of  statistical 
significance  is  available  only  for  the  F-statistic,  where  a  rigorous  P-value  is  available  based 
on  asymptotic  approximations  (see  Shumway,  1970,  1971,  1983,  2000).  The  cross  correlation 
is  based  on  the  lag  corresponding  to  the  local  maximum,  a  function  of  a  highly  variable 
measurement  that  ignores  the  fact  that  one  must  use  cross  correlations  involving  common 
sensors.  The  Capon  estimator  has  a  distribution  (see  Capon  and  Goodman,  1979)  that 
depends  on  the  noise  spectrum,  a  nuisance  parameter  that  would  have  to  be  estimated  from 
noise  preceding  a  potential  signal.  The  MUSIC  estimator  depends  on  deriving  distribution 
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theory  for  a  statistic  depending  on  the  eigen  vectors  of  the  spectral  niatrix  and  will  be 
unnecessarily  complex  (see  Stoica  and  Nehorai,  1989),  given  the  erratic  performance  that 
the  method  exhibits  when  applied  to  determining  the  exact  number  of  signals.  Nevertheless, 
the  array  version  of  this  statistic  has  not  been  applied  in  a  seismic  context,  and  we  retain  it 
as  an  option. 

4.  Multiple  Signal  Detection  and  Estimation 

It  is  obvious  from  the  preceding  section  that  the  conventional  methods  are  unreliable  for 
detecting  the  correct  number  of  signals  and  for  estimating  their  velocities  and  azimuths. 
In  this  section,  we  develop  a  sequential  extension  of  the  F-statistic  for  multiple  signals 
that  provides  effective  methodology  for  identifying  the  correct  number  of  signals  and  their 
estimated  velocities  and  azimuths.  We  also  provide  estimators  for  the  uncertainty  using  a 
multiple  signal  bootstrap  procedure.  Finally,  we  show  how  to  develop  deconvolved  waveforms 
for  each  of  the  component  signals. 

Solutions  to  the  above  problems  would  enable  online  detection  procedures  that  would  identify 
single  signals  and  would  also  identify  multiple  signals  or  contaminating  propagating  noises. 
Such  mixtures  contribute  to  errors  in  arrival  time  estimation  that  is  essential  for  location 
purposes.  Uncertainty  estimators  for  velocity  and  azimuth  would  help  to  identify  travel 
time  uncertainties  in  un-calibrated  areas.  Procedures  for  estimating  velocities  of  P  and  S 
components  could  be  helpful  in  constructing  travel  time  predictions  for  separate  phases  for 
regional  data.  Deconvolutions  could  be  used  to  determine  better  amplitudes  and  spectral 
ratios  for  discrimination  purposes.  It  should  be  noted  that  work  on  multiple  signal  models 
has  been  common  in  early  seismic  research  (see  Shumway  and  Dean,  1968,  Shumway,  1970, 
1971,  Blandford  et  al  1973,  Shumway,  1983). 

4.1  Maximum  Likelihood  Estimation 

The  proposed  solutions  to  the  above  problems  rest  squarely  on  formulating  the  multiple 
signal  model  as  a  nonlinear  regression  problem  in  the  frequency  domain.  We  formulate  a 
model  for  a  time  series  observed  at  sensor  j  =  N  sA,  time  t  =  1, . . . ,  n  as  the  delayed 
sum  of  M  signals  and  noise,  of  the  form  (1).  In  particular,  we  wish  to  have  an  unambiguous 
sequential  procedure  for  estimating  the  number  of  components  M  and  for  estimating  the 
slowness  matrix  0  for  a  given  value  of  M.  In  this  analysis,  the  signals  are  assumed  to  be 
deterministic  unknown  functions,  and  the  noises  axe  assumed  to  be  stationarily  correlated 
over  time  but  independent  from  channel  to  chaimel.  As  usual,  we  work  in  the  frequency 
domain  by  vectorizing  yj{t)  and  taking  Discrete  Fourier  Transforms  to  obtain  the  model 
(2).  In  general,  we  evaluate  all  statistics  over  L  =  DT  frequencies,  where  B  denotes  the 
bandwidth  of  interest  in  Hz  and  T  denotes  the  sample  length  in  seconds.  The  AT  x  1  vector 
Yt  is  expressed  in  terms  of  the  unknown  signal  transform  Si  and  the  NxM  slowness  matrix 
Zi{Q)  as  in  (2)  and  (3).  The  unknown  slowness  parameters  are  collected  in  the  2  x  M  matrix 
©  =  (01, 02  •  •  • ,  Gm)-  Regression  theory  (see  Shumway  and  Stoffer,  2000,  Section  5.4)  applied 
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for  a  given  slowness  matrix  0  implies  that  we  may  estimate  the  rmknown  signal  by 


5,(0)  -  [^;(0)^,(0)]  "zmYe. 


Generally,  5,(0)  is  concentrated  out  of  the  likelihood,  and  we  determine  0  as  the  minimizer 
of 

SSE{Q)  =  X:  WiYi  -  2,(0)5,(0))||^  (19) 

Searching  for  a  minimum  in  the  squared  error  (19)  over  0  or,  equivalently,  maximizing  the 
F-statistic  for  testing  5,(0)  =  0  leads  to 


F(0)  =  c 


Pb{Q) 

ssE{ey 


where 


Pb(0)  -  X]y,*z,(©)b;(©)2,(©)]'''z;(0)r, 

r=i 


is  the  generalized  beam-power  and  c  —  df2/dfi  is  the  constant  required  to  convert  (20)  to  an 
F-statLstic  with  dfi  =  2(5  -f  M)  and  dfa  =  2[5(Ar  —  M)  —  1]  degrees  of  freedom.  We  call  (21) 
the  generalized  beam  power  because  it  reduces  to  the  ordinary  beam  power  when  M  =  1. 

4.2  Sequential  Analysis  Using  Partial  F-statistics 

In  order  to  develop  a  sequential  F-test  for  adding  in  the  signal  components  81,82,  ■  ■. ,  8l, 
we  appeal  to  the  likelihood  ratio  test  of  5,  =  (5,i,0)  against  the  alternative  5,  =  (5, 1,5,2) 
where  the  signal  vector  is  partitioned  into  Afi  and  M2  components  {M  —  Mi  +  M2)  respec¬ 
tively.  The  likelihood  ratio  criterion  yields 


F==c 


555(0i,O)-55g(0i,02) 

555(01,02) 


where  c  ==  df2/dfi  and  ©1  denotes  the  estimator  for  the  slowness  matrix  ©i  under  the  reduced 
model.  The  test  statistic  compares  the  mean  squared  error  under  the  hypothesis  that  the 
added  signal  is  absent  under  the  h)q)othesis  that  the  added  signal  is  present  (numerator)  with 
the  denominator  which  measures  the  mean  squared  error  under  the  full  model.  Under  the  null 
hypothesis,  F  has  an  F  distribution  with  with  dfi  —  2M2{L  +  1)  and  d/2  =  2[L{N~M)  —  1)] 
degrees  of  freedom.  To  apply  the  result,  we  add  the  potential  signals  in  sequentially,  each 
time  taking  Mi  as  the  number  currently  in  the  model  and  M2  =  1.  It  is  intuitively  pleasing 
that  the  numerator  is  essentially  the  difference  between  the  generalized  beam  power  under 
the  reduced  model  and  the  full  model.  A  sequential  search  proceeds  by  performing  the  partial 
F-test  (22)  with  M  taking  values  1,2,3,...  and  stopping  when  F  is  no  longer  statistically 
significant. 
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Note  that  a  different  sequential  procedure,  based  on  residuals  from  lower  order  signal  models, 
was  proposed  in  Smart  (1972)  who  searched  for  a  second  signal  by  first  stripping  the  first 
signal  away,  i.e.,  by  computing 

V£i  —  Yt  —  Zii{QiSn 

and  then  selecting  ©2  as  the  value  minimizing 

Y,\\Vn-Z£{Q2)Sn\?. 

t=\ 


where  Sn  denotes  the  solution  of  (18)  under  the  reduced  model  and  5^2(©2)  denotes  the 
solution  of  (18)  with  the  residuals  V t\  as  observations. 

An  additional  statistic  of  interest  in  applications  is  the  proportion  of  power  accounted  for 
by  any  given  multiple  signal  model,  given  as 


R^{e)  = 


Pb{Q) 


(23) 


which  is  immediately  recognized  as  the  ratio  of  the  generalized  beam  power  to  the  total 
power.  Smart  has  suggested  using  the  proportion  of  power  accounted  for  as  a  measure  of 
how  well  the  model  fits. 


4.3  Model  Selection 


An  alternate  approach  can  be  taken  in  terms  of  a  model  selection  criterion  in  the  AIC  family. 
We  take  here  a  modification  to  the  corrected  AIC,  developed  especially  for  regression  models 
by  Hurvich  and  Tsai  (1989).  In  the  case  of  the  signal  model  (2),  we  develop  a  real  regression 
model  isomorphic  to  (2)  and  then  expand  in  a  Taylor’s  serias  about  the  2M  x  1  vector  9. 
Taking  the  Hurvich-Tsai  approach  essentially  doubles  the  numbers  appearing  in  the  usual 
form.  We  obtained 


AICc{M) 


,  SSE{Q) 
log  -  -  -  ^  + 


2NL  +  2M{L  +  1) 


NL  2NL--2M{L  +  l)-2' 


(24) 


The  value  is  computed  for  signal  models  with  respectively  M  ~  1,2,3,...  and  the  value  M 
that  minimizas  AIC{M)  determines  the  number  of  signals  chosen  for  the  model. 


As  a  confirming  check,  we  also  compute  the  values  of  the  model  selection  criterion  (17)  as  a 
function  of  m  and  look  for  a  confirming  minimum.  We  illustrate  the  approach  for  the  long 
and  short  period  events  given  in  Figures  1  and  3  and  show  how  it  leads  to  resolving  the 
mixtures  of  signals  for  these  two  data  sets. 


4.4  Confidence  Intervals 


First,  note  that  we  have  ended  with  a  multiple  signal  model  and  need  the  variances  and 
covariances  of  the  slowness  matrix  0  =  •  •  •  i^m)-  For  the  vector  case,  it  would  be 
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possible  to  express  the  non-linear  regression  model  (2)  as  a  first-order  Taylor’s  series  ex¬ 
pansion  in  ©  and  81,82, ...  ,8l.  Then,  the  Gauss-Newton  covariance  matrix  derived  firom 
the  nonlinear  regression  analysis  prevails.  The  variance  covariance  matrix  of  the  estimated 
velocity  and  azimuth  parameters  could  then  be  computed  by  the  delta  method.  This  proce¬ 
dure  seems  rather  involved,  and  we  adopt  a  simpler  approach  that  depends  on  a  frequency 
domain  version  of  the  bootstrap  (Paparoditis  and  Politis,  1999,  Shumway  and  Stoffer,  2000, 
p.  245).  Note  that  the  residuals 

Vi=Ye-Zeie)8S)  (25) 

should  be  approximately  independent  with  mean  zero  and  variance  fv(ui)  under  the  model. 
In  order  to  obtain  one  bootstrap  sample,  take  a  random  sample  of  the  residuals  above,  with 
replacement.  Suppose  this  random  sample  is  denoted  by  Vf  \  . . .,  For  this  initial 

bootstrap  sample,  construct  a  pseudo-sample  of  observations 

=  Zt{Q)8S)  +  Vf^  (26) 

and  apply  the  estimation  procedure,  obtaining  an  estimated  0^^^  for  this  first  bootstrap 
sample.  Then,  repeat  the  above  sampling  procedure,  obtaining  a  large  number,  say  500, 
estimated  values  ©^^1, . . . ,  Finally,  convert  these  values  to  velocity  and  azimuth 

and  use  the  two  sampling  distributions  to  obtain  variances  and  lower  and  upper  confidence 
intervals  for  velocity  and  azimuth. 

4.5  Deconvolution  of  Multiple  Signals 

It  should  be  noted  that  the  waveforms  of  the  signals  recovered  from  the  final  model  may  show 
features  that  are  not  available  in  the  simple  beam.  The  signals  may  be  recovered  using  the 
inverse  finite  Fourier  transform  of  the  frequency  domain  version  (14)  which  will  be  available 
over  a  hmited  band  spanned  by  the  vector  signal.  In  order  to  estimate  the  vector  of  time 
functions  s(t)  —  . . .  ,SM{t)y,  we  expand  the  the  fi-equency  range  to  (—1/2, 1/2)  or 

(0, 1)  by  taking  8{oJe}  =  Sf,  over  the  band  and  8{u}£)  —  0  for  uii  —  ijn,  ^  =  0, 1, . . . ,  n/2  and 
by  completing  to  frequencies  between  1/2  and  1  using  5(a;f+„/2)  =  conj  8{^t...ni2),.  Taking 
the  inverse  Fourier  transform 

n— 1 

s{t)  —  S{u}t)exp{2mujei}  (27) 

e=o 

and  shifting  the  coefficients  gives  the  deconvolved  series.  Note  that  increasing  the  original  n 
to  2n  will  eliminate  the  end  effects  of  the  discrete  Fourier  transform. 

5.  Multiple  Signal  Analysis  of  Recorded  Events 

The  value  of  multiple  signal  techniques  developed  in  Section  4. 1-4.3  over  and  above  conven¬ 
tional  multiple  signal  techniques  such  as  MUSIC  or  single  signal  techniques  like  the  F-statistic 
(semblence)  or  the  Capon  estimator  ultimately  rests  on  their  success  in  analyzing  real  data. 
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Motivation  for  the  multiple  signal  analysis  is  provided  initially  by  analyzing  a  teleseismic 
recording  known  to  contain  two  earthquakes  where  it  is  shown  that  conventional  techniques 
fail  and  that  a  combination  of  the  partial  F-statLstic  and  model  selection  separate  the  record¬ 
ing  into  two  earthquakes  and  a  noise  source,  known  to  be  generated  by  a  Pacific  storm.  Ver¬ 
ification  of  the  methodology  is  provided  by  using  a  contrived  mixture  of  signals  with  known 
velocities  and  azimuths.  A  third  example  involves  analyzing  a  noisy  China  event  where  the 
elimination  of  propagating  noiseg  provides  a  visual  verification  of  an  obscured  depth  phase. 

While  the  examples  make  it  clear  that  conventional  methods  work  poorly,  that  is  not  the 
only  motivation  for  using  the  techniques  derived  in  Section  4.  We  mention  the  following 
specific  advances  for  the  new  methodology: 

A.  Nonlinear  optimization  in  the  sequential  procedure  provides  computationally  more 
accurate  velocities  and  azimuths  than  are  obtained  by  grid  search  methods  typically 
used  by  conventional  methods  in  Section  3.  All  velocities  and  azimuths  are  estimated 
by  maximum  likelihood  and  hence,  will  be  efficient. 

B.  The  sequential  procedure  using  partial  F-statistics  and  AICc  shows  the  results  for  all 
possible  multiple  signal  configurations  and  settles  on  the  one  that  minimizes  AICc 

C.  Uncertainties  for  the  best  model  in  2.  are  computed,  giving  standard  deviations  and 
confidence  intervals  for  velocity  and  azimuths. 

D.  The  separate  waveforms  of  the  component  signals  are  estimated  using  least  squares 
deconvolution  in  the  frequency  domain. 

We  should  make  it  clear  from  the  onset  that  the  same  systematic  procedure  for  identifying 
the  components  of  the  possible  mixture  and  estimating  the  resulting  velocities  and  azimuths 
of  the  components  will  be  followed  in  each  of  the  examples  in  the  following  sections.  We 
summarize  the  approach  for  each  example  as: 

1.  An  initial  time~frequency  spectral  analysis  is  used  to  determine  an  appropriate 
bandwidth  for  the  conventional  and  sequential  analysis.  This  analysis  can  also  be  used 
to  focus  on  particular  arrival  phases  and  to  isolate  the  appropriate  time  interval  and 
scale  for  the  sequential  analysis. 

2.  Conventional  slowness  analysis  plots  for  (a)  beam  power,  (b)  Capon,  (c)  F- 
statistic  and  (d)  MUSIC  estimators  are  used  to  suggest  regions  in  slowness  space 
where  there  might  be  signals.  This  will  often  help  in  setting  initial  slowness  values  for 
the  multiple  signal  optimization  sequence. 

3.  An  initial  maximum  for  the  number  of  signals  (usually  no  more  than  four)  and  their 
initial  slowness  coordinates  is  used  to  start  a  sequential  analysis,  with  the  final  momber 
of  signals  chosen  by  minimizing  the  fit  Corrected  Information  Theory  Criterion 
AICc-  The  proportion  of  variability  accounted  for  and  the  partial  F-statistic 
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for  adding  a  signal  are  shown  at  each  step.  The  final  configuration  is  determined  by  a 
subjective  weighting  of  these  criteria. 

4.  The  frequency  domain  bootstrap  is  used  to  develop  95%  confidence  intervals  for  the 
velocities  and  azimuths  in  the  optimal  configuration. 

5.  Deconvolution  is  used  to  produce  estimated  waveforms  for  the  component  signals. 

5.1  A  Teleseismic  Recording  Containing  Two  Earthquakes  and  a  Noise  Source 

We  begin  by  considering  the  data  that  precipitated  interest  in  multiple  signal  models.  This 
data,  shown  in  Figure  1,  contains  a  verifiable  mixture  of  two  earthquakes,  one  from  the  south 
of  Africa  and  another  from  the  Philippines  observed  at  the  Korean  Seismic  Array  (KSAR). 
Figure  2  repeats  the  last  channel  of  this  mixture  and  shows  its  time-varying  spectrum. 

The  time-frequency  spectrum  suggests  a  band  (.02-.08  Hz)  for  analysis  and  we  adopted  this 
range  for  the  slowness  statistics  shown  in  Figure  3.  The  time-frequency  spectrum  shows 
two  main  arrivals  with  2-3  possible  secondary  arrivals.  Figure  3  shows  the  slowness  plots 
for  methods  (a)-(d)  and  we  note  that  the  beam  power  and  F-statistic  show  only  a  single 
arrival  at  203  degrees.  All  methods  show  initial  arrivals  in  this  range,  with  the  Capon 
and  MUSIC  (3  signals  assumed)  estimators  showing  bulges  in  the  150-160  degree  range.  The 
correct  azimuths  are  226  degrees  and  198  degrees  which  are  both  different  than  the  estimated 
primary  azimuth  obtained  by  all  methods.  The  results  are  smnmarized  in  Table  1  below. 

Table  1;  Conventional  Estimates  for  Long  Period  Event 


Estimates 

F-stat 

Correlation 

MUSIC 

Capon 

Azimuth  1 

203 

204 

204 

207 

Velocity- 1 

3.5 

3.7 

3.7 

4.1 

Azimuth-2 

158 

149 

Velocity-2 

5.3 

4.5 

We  complete  Step  3.  by  setting  an  initial  guess  for  the  maximum  signals  at  4,  with  initial 
guesses  for  the  nonlinear  optimization  focusing  first  on  the  third  quadrant  (-.1  -.1),  with  a 
possible  secondary  in  (.1,  -.1).  The  other  two  start  slowness  values  were  (.1,  .1)  and  (-.1,  .1). 
Table  2  below  shows  the  results  of  the  model  selection  procedure  and  the  results  indicate 
three  signals  by  AICc  and  possibly  four  signals  by  the  partial  F-test.  Values  around  2  are 
marginal  for  this  particular  F  and  we  settle  on  three  signals. 

Table  2:  Model  Selection  for  Long  Period  Event 


Model 

R2 

WSBk 

Partial  F-stat 

1-Signal 

.86 

-1.25 

35.1 

2-Signal 

.95 

-1.87 

3-Signal 

.98 

-2.21 

4-Signal 

.99 

-1.57 
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Figure  2:  Teleseismic  signal  of  Figure  1  at  Korean  Seismic  Array  with 
time- varying  spectrum. 
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Figure  3:  Slowness  plots  for  the  teleseismic  signal  in  Figure  1.  The  circle 
denotes  a  velocity  of  3km/sec. 
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The  final  results,  including  the  uncertainty  estimates  using  the  bootstrap  as  in  Step  4  of  the 
suggested  procedure  leads,  are  summarized  in  Table  3  below.  We  note  that  the  estimators 
are  still  several  degrees  off  the  known  azimuths  but  that  the  confidence  intervals  include  the 
true  values. 

Table  3:  Summary  Estimators  and  Uncertainties  for  Long  period  Event 


Parameter 

Known  Value 

Estimate(sd) 

95%  Confidence 

Azimuth  1 

198 

iQm 

197-203 

Azimuth  2 

- 

125-134 

Azimuth  3 

226 

216-229 

The  final  step  involves  deconvolving  the  signals  and  the  separate  components  are  shown  in 
Figure  4.  The  signal  at  223  degrees  looks  reasonable  for  the  primary  signal,  as  does  the 
second  signal  at  200  degrees.  The  noise  at  130  degrees  seems  large  and  there  must  be  some 
cancelation  for  the  three  to  add  to  the  observed  series.  For  this  reason,  it  seems  essential 
to  construct  a  mixture  where  the  component  signals  are  known  and  we  do  this  in  the  next 
section. 


Figure  4:  Deconvolution  of  two  identified  signals  at  known  azimuths  of  226 
and  198  degrees  and  a  third  unidentified  signal  or  noise  source 
at  135  degrees  at  a  long  period  array. 
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5.2  Two  Regional  Events  and  a  Contrived  Mixture 

In  order  to  verify  that  the  procedures  described  in  Sections  4. 1-4.5  perform  as  claimed,  we 
took  regional  signals  recorded  at  the  Korean  seismic  array  from  38  degrees  and  196  degrees 
respectively  and  first  analyzed  them  separately  using  Steps  1-4.  Figure  5  shows  the  time- 
frequency  spectrum  of  the  regional  signal  and  suggests  a  frequency  range  of  .3-3.3  Hz  for 
analysis.  The  first  10  seconds  seems  to  contain  a  mixture  of  three  different  frequencies, 
centered  at  about  .75,1.5  and  2.8  Hz  respectively.  The  signal  is  quite  impulsive.  Figure  6 
shows  a  strong  arrival  at  38  degrees  (F-statistic  and  cross  correlation)  or  39  degrees  (Capon 
and  MUSIC).  The  results  are  summamrized  in  Table  4  below. 

Table  4:  Conventional  Estimates  for  Regional  Signal  from  38  Degrees 


Estimates 

F-stat 

Correlation 

MUSIC 

Capon 

Azimuth 

38 

39 

39 

Velocity 

11.9 

15.6 

15.6 

Step  3  can  be  performed  as  a  check  against  the  possibility  of  additional  signals.  Starting 
with  the  assumption  that  there  are  a  maximum  of  two  signals  present,  we  chose  (.1,  .1)  for 
the  first  slowness  and  (-.1,  -.1)  for  the  start  point  in  a  search  for  a  potential  additional  signal, 
we  obtained  the  results  as  summarized  in  Table  5  below.  We  note  that  AICc  is  minimized 
for  the  single  signal  model  and  that  the  partial  F-statistic  of  1.7  is  not  large  enough  to  add 
the  second  signal.  Note  also  that  the  percentage  of  variation  accounted  for  (57%)  is  much 
smaller  for  this  noisier  regional  data. 

Table  5;  Model  Selection  for  38  Degree  Regional  Signal 


Model 

R2 

AICc 

Partial  F-stat 

1-Signal 

.57 

-2.96 

22.8 

2-Signal 

.61 

-2.93 

1.7 

The  single  signal  model  then  gives  estimated  azimuths  of  38.3(.3)  degrees  with  the  95% 
confidence  interval  (37.9-39.0)  via  the  bootstrap.  The  estimated  velocity  was  11.9(.l)  km/sec 
with  95%  confidence  interval  (11.8-12.1)  km/sec.  Using  the  slowness  vector  (.0519,  .0657) 
corresponding  to  the  these  velocities  and  azimuths  leads  to  the  deconvolution  estimator 
shown  in  Figure  9. 

As  a  second  regional  event,  consider  another  regional  event,  shown  in  Figure  7,  recorded 
at  the  same  array.  This  particular  event  generated  a  longer  lower  frequency  signal,  as  can 
be  seen  from  the  single  channel  and  time-frequency  spectrum  shown  in  Figure  7.  Figure  8 
shows  the  slowness  plots  for  this  second  event  and  it  can  be  seen  that  all  methods  indicate 
a  single  peak  for  velocity  and  azimuth,  with  the  summary  r&sults  as  given  in  Table  6. 

We  note  here  the  beginnings  of  a  persistent  tendency  of  the  Capon  and  MUSIC  estimators 
to  produce  higher  than  normal  velocities.  As  can  be  seen  in  (14)  and  (15),  both  of  these 
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Figure  5:  Regional  signal  (38  degrees)  at  Korean  Seismic  Array  with  time- 
varying  spectrum. 
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Figure  6:  Slowness  plots  for  the  regional  signal  from  38  degrees  at  Korean 
Seismic  Array.  The  circle  denotes  a  velocity  of  8km/sec. 
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Figure  7:  Regional  signal  (196  degrees)  at  Korean  Seismic  Array  with 
time- varying  spectrum. 
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Figure  8:  Slowness  plots  for  the  regional  signal  from  196  degrees, 
circle  denotes  a  velocity  of  8km/sec. 
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estimators  involve  the  spectral  matrix  (13)  which  must  be  estimated  from  a  band  of  of  L  > 
frequencies  where  N  denotes  the  number  of  sensors.  The  matching  vector  zi{0)  is  chosen  at 
the  center  £  corresponding  to  the  average  frequency  over  the  interval.  It  might  be  better  to 
modify  both  statistics  by  averaging  over  the  single  frequency  matches  but  we  have  not  tried 
this  approach.  The  azimuths  for  all  methods  match  quite  nicely. 

Table  6;  Conventional  Estimates  for  Regional  Signal  from  196  degrees 


Estimates 

F-stat 

Correlation 

MUSIC 

Capon 

Azimuth 

196 

196 

198 

198 

Velocity 

9.4 

9.2 

15.8 

15.8 

Table  7;  Model  Selection  for  196  Degree  Regional  Signal 


Model 

AICc 

Partial  F-stat 

1-Signal 

.59 

-1.61 

24.3 

2-Signal 

.65 

-1.62 

7.5 

3-Signal 

.68 

-1.56 

1.4 

The  model  selection  results  for  this  particular  regional  signal  indicate  that  a  second  contar- 
minating  process  may  be  present.  The  model  selection  statistic  AICc  is  minimized  for  the 
two-signa.l  model  and  the  F-statistic  for  this  signal  may  be  too  large  to  ignore.  Our  tenta¬ 
tive  conclusion  for  this  event  is  that  there  is  a  primary  signal  at  196(.6)  degrees  with  95% 
confidence  interval  (195-197)  degrees  and  a  secondary  signal  at  56(1.1)  degrees  with  95% 
confidence  interval  (55-59)  degreas.  The  velocity  of  this  second  signal  is  slower  at  3.1  (.05) 
km/sec.  However,  the  joint  deconvolution  of  the  two  components  using  the  slowness  vector 
resulting  from  the  two-signal  model  is  shown  in  Figure  10.  Clearly,  the  results  show  that  the 
contaminating  signal  has  low  amplitude  and  we  note  parenthetically  that  it  probably  will 
not  exert  a  significant  effect  on  the  contrived  mixture. 

As  a  test  of  the  estimation  procedure,  the  two  previous  arrays  were  mixed  by  adding  the  two 
events  together  on  common  channels  with  equal  amplitudes.  Figure  11  shows  the  contrived 
mixture  along  with  its  time-frequency  spectrum.  Note  first  that  the  time-frequency  spectnim 
suggest  the  frequency  range  .2-3.3  Hz  and  this  is  what  was  used  as  input.  The  slowness  plots 
all  show  both  signals  as  being  present  at  the  correct  azimuths.  Note  that  the  Capon  and 
MUSIC  again  over-estimate  the  velocities  if  the  average  frequency  of  1.75  Hz  is  chosen. 
Choosing  .9  Hz  on  the  basis  of  the  time-frequency  spectrum  leads  to  the  more  reasonable 
estimates  shown  in  the  plot  and  in  Table  8. 
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Table  8:  Conventional  Estimates  for  Contrived  Regional  Mixture 


Estimates 

F-stat 

Correlation 

MUSIC 

Capon 

Azimuth- 1 

196 

196 

196 

196 

Velocity- 1 

9.5 

9.4 

13.7 

13.7 

Azimuth-2 

38 

38 

38 

Velocity-2 

9.9 

7.5 

7.5 

It  was  noted  in  Figure  12  that  the  slowness  plots  for  the  Capon  and  MUSIC  (with  two  signals 
assumed)  detectors  showed  the  second  weaker  signal  to  be  present  at  fairly  clearly  delBned 
azimuths.  While  this  simple  case  is  one  where  a  weaker  but  well  separated  signal  is  detected  by 
these  two  methods,  tests  of  significance  and  confidence  intervals  are  not  provided.  The  single 
signal  version  of  the  F-statistic  shows  only  a  weak  indication  of  a  second  signal  so  that  it  is  clear 
that  more  must  be  done. 

Hence,  we  apply  the  partial  F-sta.tistic  and  AICq,  fitting  a  sequence  of  models  with  M  — 
1, 2, 3  signals.  The  results,  summarized  in  Table  9  below,  show  the  correct  number  of  signals 
by  either  the  model  selection  or  partial  F-statistic  criteria.  Therefore,  the  two-signal  model 
can  be  used  to  develop  the  final  estimators  and  confidence  intervals  for  velocity  and  azimuth 
which  will  produce  the  deconvolution  of  the  two  signals. 

Table  9:  Model  Selection  for  Contrived  Regional  Mixture 


Model 

R2 

AICc 

Partial  F-stat 

1-Signal 

.46 

-1.90 

14.5 

2-Signal 

.60 

-2.05 

5.3 

3-Signal 

.64 

-2.00 

1.6 

Table  10  summarizes  the  final  estimators  and  their  uncertainties  for  the  optimum  two-signal 
model.  Note  that  both  azimuth  estimators  are  close  to  those  obtained  for  the  separate  event 
analyses  so  this  simple  example  verifies  that  the  methodology  used  for  the  long  period  event 
is  reasonable. 

Table  10:  Summary  Estimators  and  Uncertainties  for  Contrived  Mixture 


Parameter 

Known  Value 

Estimate(sd) 

95%  Confidence 

Azimuth  1 

196 

195.5(.62) 

194.1-196.5 

Velocity  1 

9.5(.ll) 

9.2-9.7 

Azimuth  2 

38 

38.4(.88) 

37.8-41.2 

Velocity  2 

12.0(.14) 

11.8-12.3 

Finally,  Figure  13  shows  the  result  of  the  least  squares  deconvolution  of  the  two  signals. 
Comparing  these  deconvolved  signals  with  the  single-event  deconvolutions  in  Figures  9  and 
10  shows  excellent  agreement. 
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Figure  13:  Deconvolution  of  mixture  (38,  196  degrees),  degrees. 


5.3  A  Regional  Event  from  China  With  Depth  Phase  Obscured  by  Noise 

In  this  example,  we  look  further  at  the  event  analyzed  by  Stroujkova  and  Rieter  (2006)  in 
their  preliminary  report.  The  event,  from  NE  China,  had  magnitude  of  4.7,  with  published 
depths  between  0  and  15  kilometers.  Because  of  the  high  noise  level,  as  can  be  noted  from 
Figure  14,  there  will  be  difficulty  in  establishing  whether  or  not  there  is  a  depth  phase.  We 
show  in  this  example  how  a  multiple  signal  approach  eliminates  enough  of  the  noise  so  that 
the  deconvolved  signal  shows  a  clear  depth  phase. 

The  time-frequency  analysis  in  Figure  14  covers  a  relatively  long  time  interval.  The  broader 
time  interval  is  primarily  to  emphasize  the  multiple  arrivals  and  relatively  complex  structure 
of  the  single  channel  time-frequency  spectra. 

Figure  15  shows  the  slowness  analysis  for  this  event  and  we  note  that  the  conventional 
estimators  show  a  primary  arrival  in  the  neighborhood  of  284  degree  except  for  the  correlation 
analysis  which  is  off  by  about  40  degrees.  Furthermore,  the  Capon  and  MUSIC  estimators 
show  a  secondary  peak  in  the  150-160  degree  range. 
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Figure  14:  The  1991  China  event  and  its  time  varying  spectrum. 
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Figure  15:  Slowness  plots  for  the  1991 
velocity  of  8km/sec. 
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Table  11:  Conventional  Estimates  for  Regional  Event  From  China 


Estimates 

F-stat 

Correlation 

MUSIC 

Capon 

Azimuth-1 

287 

247 

284 

286 

Velocity- 1 

8.1 

24.9 

12.1 

13.7 

Azimuth-2 

153 

159 

Velocity-2 

13.5 

14.1 

The  results  of  a  combined  multiple  signal  analysis  using  the  partial  F-statistics  and  the  model 
selection  criterion  AICc  are  shown  in  Table  12.  Note  that  two  signals  each  give  statistically 
significant  added  power,  as  measured  by  the  partial  F-statistics  and  that  the  third  one  does 
not  add  substantially  to  the  power.  Noise  is  still  quite  high,  with  only  50%  of  the  power 
accounted  for  by  the  two  signals,  accounted  for  by  the  two  main  signals.  The  estimated 
velocities  and  azimuths  of  the  two  signals  are  287  degrees  and  8.2  km/sec  for  the  first  and 
151  degrees  and  3.4  km/sec  for  the  second  signal.  The  estimated  velocity  of  the  second  .signal 
differs  from  that  estimated  by  MUSIC  but  some  runs  at  different  frequency  centers  yielded 
estimators  in  that  range  for  the  other  estimators  as  well.  The  estimators  are  sununarized  in 
Table  13. 


Table  12:  Model  Selection  for  China  Event 


Model 

R2 

AICc 

Partial  F-stat 

1-Signal 

.34 

-.40 

8.1 

2-Signal 

.46 

-.43 

3.4 

3-Signal 

.52 

-.41 

1.6 

Table  13:  Summary  Estimators  and  Uncertainties  for  Noisy  China  Event 


Parameter 

Estimate(sd) 

95%  Confidence 

Azimuth  1 
Velocity  1 
Azimuth  2 
Velocity  2 

150.8(1.3) 

3.4(.08) 

286.8(1.1) 

8.2(.ll) 

148.2-153.8 

3.3-3.6 

284.5-288.4 

8.0-8.5 

The  deconvolution  corresponding  to  the  two  components  noted  in  the  previous  section,  shown 
in  Figure  16,  gives  the  estimated  signal  and  noise  components  in  the  regional  data.  The 
estimated  noise  component  from  148  degrees  shows  the  character  of  the  low  frequency  noise 
before  the  signal  enters.  The  estimated  signal  from  286  degrees  shows  the  depth  phase  pP 
much  more  clearly  than  does  the  single  channel  mixture.  In  this  case,  one  obtains  from  the 
plot  a  delay  of  about  4.6  seconds  which  is  comparable  to  that  obtained  by  Stroujkova  and 
Reiter  (2006).  We  also  note  the  noise  reduction  capabilities  of  the  two-signal  beam. 
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Figure  16:  Deconvolution  of  Chinese  event  into  a  noise  source  and  a  287 
degree  signal. 


6.  Software  Documentation  and  Data  Files 

6.1  General  Structure 

The  software  delivered  with  this  final  report  is  structiired  as  a  main  driving  program  with 
function  subroutines  for  computing  and  plotting  all  quantities  given  in  Sections  3,  4  and  5  of 
this  report.  The  one  exception  is  the  computation  and  plotting  for  the  original  series  and  its 
time^-frequency  spectrum  which  we  assume  will  be  available  to  any  user.  The  main  program, 
as  indicated  below  calls  the  five  functions  that  are  needed  to  perform  the  multiple  signal 
a.nalysis. 

The  function  beam.m  plots  four  conventional  estimators  (beam,  Capon,  F-statistic  and  MU¬ 
SIC)  in  slowness  space  and  gives  the  estimated  velocities  and  azimuths.  A  separate  function 
corr_sl.m,  gives  the  broad-band  correlation  estimators.  For  the  contrived  mixture  of  signals 
shown  in  the  text,  the  two  functions  produce  Figure  12  and  the  entries  necessary  for  Table 
8. 

The  function  vmiLsig.m  develops  the  multiple  signal  analysis  described  in  Section  4.  This 
gives  the  sequence  of  partial  F-statistics  and  values  of  AICc  necessary  to  settle  on  on  op¬ 
timum  number  of  signals.  The  function  identifies  the  optimum  number  of  signals  and  gives 
the  corrasponding  slownesses  for  input  to  the  deconvolution.  The  model  selection  criteria 
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are  computed  for  inclusion  in  Table  9.  The  estimators  for  velocities  and  azimuths  corre¬ 
sponding  to  the  optimum  model  are  given  and  the  option  to  compute  standard  errors  and 
95%  confidence  intervals  is  included  in  the  function  confJnt.m.  This  would  complete  the 
entries  for  Table  10.  Finally,  the  option  for  deconvolving  the  component  signals  is  given  by 
decon.m.  We  summarize  the  structure  below  where  each  main  function  is  identified  with  the 
appropriate  section  of  the  report: 


Main  Program:  Msignal 

1.  function6eam.m(3.1-3.3) 

function^rob.m 
function  m 

2.  function  corr.sl'm(3A) 

3.  function  muLsig.m{4.2-4.3) 

function  fprob.m 
function  v^az.m 
function  SSEL.m 

4.  function  decon.m{4.5) 

5.  function  confJnt.m(4.4) 

function  m 

function  SSEL.m 


Note  that  there  are  three  additional  functions  that  are  used  by  certain  some  of  the  five  main 
functions.  The  function  fprob.  m  computes  the  P-value  for  the  F-statistics  obtained  at  various 
stages.  The  function  v.az.m  computes  the  velocity  and  azimuth  corresponding  to  a  given 
slowness.  The  squared  error  objective  function  (19)  in  computed  in  the  subroutine  SSEL.m. 

6.2  Data 

We  summarize  the  data  files  provided  in  the  report  in  Table  14  along  with  the  sections  where 
an  analysis  is  provided.  More  details  for  the  inputs  are  given  in  Section  3.3 


Table  14:  Data  Provided  with  Report 


Description 

Section 

Data  File 

Array  Coordinates 

1.  Teleseismic 

5.1 

k-array 

r_k-array 

2.  Region(38dg) 

5.2 

ksl994202 

r_ksar 

3.  Region(196dg) 

5.2 

ksl99517 

rJksar 

4.  Mixture 

5.2 

ks9495mix 

r_ksar 

5.  China  Regional 

5.3 

ksl991101 

r_ksar_1991-101 
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6.3  Main  Program  and  Function  Subroutines 

The  code  for  the  main  program  is  given  below  with  the  calls  to  the  function  subroutines. 

The  inptits  are  the  data  files  and  sensor  locations,  along  with  the  sample  rate,  frequency  and 
slowness  ranges  and  an  initial  slowness  vector  for  the  largest  model  to  consider. 

In  general,  starting  the  vector  with  four  bivariate  slowness  coordinates,  i.e.  slJfull=(.sIii,  slu,  SI21 ,  .s/22,  ■■•sl 

has  been  fairly  effective.  The  single-signal  model,  M  —  1  uses  the  first  two  coordinates  as 

start  values,  the  M  =  2  signal  model  usas  the  first  two  pairs  and  the  M  =  4  model  uses  all 

four  pairs  as  start  value.  Best  results  were  obtained  by  starting  with  vectors  in  the  middle 

of  the  quadrant,  e.g.  (sn,.S2i)  =  (-I,-!)-  The  recommended  pattern  will  be  clear  for  the 

example  in  the  next  section. 

%  M_Signal.m 

%  Multiple  signal  analysis 

%  Reads  Data  and  Controls  the  Separate  Analysis  Features 

7,  Conventional  (Beam,  Capon,  F-statistic,  MUSIC,  Correlation) 

°/,  Sections  3. 1-3.4  of  Shumway  (2006) 

7,  Multiple  Signal  Analysis  (Partial  F,  AIC,  7«  Variation) 

7o  Sections  4. 1-4.3  of  Shumway  (2006) 

7o  Bootstrap  Confidence  Intervals 
7o  Section  4.4  of  Shumway  (2006) 

7.  Deconvolution 

7.  Section  4.5  of  Shumway  (2006) 

7.  Input  Information 

%  Data  File  (data):  n(points)  X  N(sensors)  rectangular  array 
7»  Sensor  Locations  (r) :  N(sensors)  X  2  Coordinates  (E,  N) 

7o  Sample  Rate  (sr) :  points/s 

7o  Frequency  Band:  (LO,  HI)  in  Hz 

7o  Slowness  Values:  Range  for  s=(s_l,s_2) 

7«  Initial  Slownesses  for  Multiple  Signal  Analysis  (sl_full) 

7»  rd=l/V_c,  V_c=constant  velocity  circle  in  slowness  plot 

id=input(’File  Number') 

M=input(' Assumed  number  of  signals  for  MUSIC  estimator') 

7»  Set  the  maximum  lag  for  the  cross  correlation  estimator 
mxlag=input(' Maximum  lag  in  pts  for  correlation  estimator  (25)') 

7«  Load  Data 
if  id==l 

load  k_axray.dat,  data=k_ array; 

7.  Sampling  rate,  frequency  range 
sr=l,L0=.02,HI=.08 
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%  Slowness  range 

s=-.5: .01: .5; ,rd=l/8; 

y.  Array  Coordinates  (km)  East-North 

sl_full=[-.l  -.1  .1  -.1  .1  .1  -.1  .1] 

load  r_k_array,  r=r_k_array; 


if  id==2 

load  ksl994202,  data-ksl994202; 

7o  Sampling  rate,  center  frequency,  bandwidth 

sr=20,L0=.3,  HI=3.3 

%  Slowness  Range 

s=-. 25:. 01:. 25;,  rd=l/8; 

y,  Array  Coordinates  (km)  East-North 

load  r_ksar,  r=r_ksar; 

sl_full=[.l  .1  -.1  -.1] 

end 


if  id==3 

load  ksl99517,  data=ksl99517; 

y.  Sampling  rate,  center  frequency,  bandwidth 

sr=20,L0=.2,HI=2.0 

y.  Slowness  Range 

s=-. 25:. 01:. 25;,  rd=l/8; 

%  Array  Coordinates  (km)  East-North 
load  r_ksar,  r=r_ksar; 
sl_full=[-.l  -.1  .1  .1  .1  -.1  ] 

end 

if  id==4 

load  ks9495mix,  data=ks9495mix; 

%  Sampling  rate,  frequency  range 

sr=20,fo=l,  L0=.2,  HI=3.3 

y.  Slowness  Range 

s=-.25: .01: .25; ,  rd=l/8: 

y.  Array  Coordinates  (km)  East-North 

load  r_kscir,  r=r_ksar; 

sl_full=[-.l  -.1  .1  .1  .1  -.1] 

end 

if  id==5 

load  ksl991101 
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sr=20 

I  L0=.3;,  HI=3.0 
L0=.2;,  HI=1.5 

data=ksl991101 (3501:4524, : ) ; 

%  data=ksl991 101 (2501 : 3524, : ) ; 
rd=l/8; 

s=[-.25: .01: .25] ; 

I  s=[-.4:.01:.4]: 
sl^full=[.l  -.1  -.1  .1  -.1  -.1  ] 
load  r_ksar_1991_101 ,  r=r_ksar_1991_101; 
end 

%  Sections  3. 1-3.3  of  Shxiniway(2006) 
beam(data,r,sr,LO,HI,s,rd,M) ; 

"/o  Section  3.4  of  Shumway  (2006) 

[az_corr ,  vel_corr]  =corr_sl (data, r , sr ,iiixlag) 

'/o  Section  4. 2-4. 3  of  Shumway  (2006) 

[M_min , az_min , vel_min , sl_min] =mul_sig (data, r , sr , LO , HI , sl_f ull) 

%  Section  4.5  of  Shumway  (2006) 

idecon=input('l  for  deconvolution,  0  otherwise')  if  idecon==l 
jdecon=input('0  for  optimum  input,  1  for  user  input') 
if  jdecon==l 

sl_min=input(' Slowness  vector') 

end 

decon(data,r,sr,L0,HI ,sl_min) 

end 

%  Section  4.4  of  Shumway  (2006) 

id_boot=input ( ' 1  if  bootstrap  Cl  0  otherwise  Bootstrap  takes  time') 
if  id_boot==l 

n_boot=input('#  of  bootstrap  reps,  500  is  suggested') 
end 

if  id_boot==l 

conf _int (data , r , sr , LO , HI , sl_min ,n_boot ) 
end 
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6.4  A  Worked  Example 

We  illustrate  a  run  of  the  program  for  the  mixture  of  signals  considered  in  Section  5.2. 
For  this  case,  we  restricted  the  start  slowness  vector  to  a  maximum  of  three  signals  to  cut 
down  on  the  printout.  For  sljfull=[-.l  -.1  .1  .1  .1  -.1],  we  looked  at  a  potential  for  M  =  3 
signals  with  the  primary  searches  in  quadrants  containing  suspected  signals.  Different  start 
configurations  can  produce  different  looking  end  results  so  the  above  start  value  should  be 
checked  against  other  permutations  of  the  two  dimensional  slowness  vectors. 

Two  plots  are  produced  which  are  not  shown.  These  are  exactly  Figures  12  and  13  with 
no  labeling  on  the  plots.  Generally,  labeling  is  data  dependent  and  can  be  done  using  the 
gtext(’laber)  instruction.  The  axes  can  also  be  changed  after  the  fact  by  statements  like 
axis([l  1024  -20000  20000]). 

»  M_signal 

File  Number4  id  =  4 

Assumed  number  of  signals  for  MUSIC 
estimator2 
M  =  2 

Starting  with  file  number  4  and  the  2-signal  assumptions  for  MUSIC 

Maximum  lag  in  pts  for  correlation  estimator  (25)25 
The  starting  value  should  cover  the  maximum  time  lag  (positive  or 
negative . 

mxlag  =  25 

sr  ==  20  fo  =  1  LO  =  0.2000  HI  =  3.3000 

Sampling  rate  and  frequency  ranges  are  typical  for  regional  event 

sl_full  = 

-0.1000  -0.1000  0.1000  0.1000  0.1000  -0.1000 
avg_frequency  =  1.7500 

The  average  of  high  and  low  cutoffs  may  not  be  the  best  for  the 
probe  vector  used  in  the  Capon  and  MUSIC  estimators 

Input  user  center  frequency. 75 

fo  =  0.7500 

R  = 
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8.3774  0.0231 

0.0231  5.5422 

F_stat  =  14.5768 

Single-signal  F-statistic  and  its  P-Value 
P_val  =  7.6690e-019 

a2_F  =  196.6992 

vel_F  =  9.5783 

az_C  =  198.4349 

vel_C  =  10.5409 

az_M  =  198.4349 

vel_M  =  10.5409 

az_corr  =  195.5305 

vel_corr  =  9.3569 

M  =  3 

Single  Signal  Model  Results 
sl_est  =  -0.0287  -0.1018 

No_Signals  =  1 

F_Azimuths  =  195 . 7325 

F_Velocities  =  9.4540 

Double  Signal  Model  Results 

sl_est  =  -0.0281  -0.1015  0.0516  0.0651 

No_Signals  =12 
F_Azimuths  =  195.4661  38.4138 
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F_Velocities  = 


9.4983  12.0364 


Triple  Signal  Model  Results 

sl_est  =  -0.0279  -0.1020  0.0514  0.0650  -0.0506 

No_Signals  =123 

F.Azimuths  =  195.2898  38.3312  189.4961 

F_Velocities  =  9.4615  12.0728  3.2607 

R.squared  =  0.4620  0.5961  0.6356 

Corrected_AIC  =  -1.8968  -2.0507  -2.0029 

Error_SS  =  281.6107  211.3944  190.7318 

Total_SS  =  523.3957 

Partial  F-statistics 

Number_signals  =  1 

F.value  =  14.5475 

Number_signals  =  2 

F_value  =  14.5475  5.2958 


Number_signals  =  3 

F.value  =  14.5475  5.2958  1.6189 

Best  Model 
M.min  =  2 

az.min  =  195.4661  38.4138 

vel.min  =9.498312.0364 

sl_min  =  -0.0281  -0;i015  0.0516  0.0651 


1  for  deconvolution,  0  otherwisel 
idecon  =  1 

0  for  optimum  input,  1  for  user  inputO 
jdecon  =  0 


-0.3025 
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2 


M  = 

nplt  =  3 

1  if  bootstrap  Cl  0  otherwise  Bootstrap  takes  timel 
id_boot  =  1 

#  of  bootstrap  reps.  500  is  suggested200 

n_boot  = 

200 


Standard  deviations  for  Azimuths  and  Velocitie 
std_Az  =  0.6476  0.9115 

std_Vel  =  0.1076  0.1475  957. 

Confidence  Intervals 
Az_025  =  193.9662  37.7199 

Az_975  =  196.4673  41.3768 

Vel_025  =  9.2418  11.7695 

Vel_975  =  9.6738  12.3464 

7.  Discussion 

This  work  was  primarily  motivated  by  the  problem  of  detecting  mixtures  of  signals  on  tele- 
seismic  and  regional  arrays.  Such  undetected  mixtures  are  shown  to  give  incorrect  velocities 
and  a.zirauths  when  treated  by  traditional  single-signal  detection  methods  based  on  cross¬ 
correlation  or  frequency  wave-number  methods.  Furthermore,  the  separation  of  interfering 
phases  and  coherent  noise  sources  should  improve  detection  statistics  and  lead  to  improve¬ 
ments  in  location  and  magnitude  estimates. 

Using  current  improved  computing  platforms  such  as  MATLAB  make  the  nonlinear  estima¬ 
tion  problems  implicit  in  multiple  signal  modeling  tractable  and  easy  to  implement.  Using 
early  work  involving  multiple  signal  estimation  (Shumway,  1970)  and  its  extensions  to  wave- 
number  methods  (Smart,  1972,  1976),  we  were  able  to  formulate  the  problem  in  a  regression 
framework  that  led  to  a  sequential  detection  approach  for  reliably  determining  the  munber  of 
signals  and  coherent  noises  present  along  with  their  estimated  velocities  and  azimuths.  We 
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have  provided  methods  for  comparing  detection  performance  using  the  F-statistic,  the  AICc 
model  selection  statistic,  and  confidence  intervals  for  the  velocities  and  azimuths  using  the 
bootstrap.  Finally,  we  show  that  the  regression  model  provides  a  method  for  deconvolving 
the  component  signals  and  exhibit  results  for  both  long  and  short  period  seismic  data. 
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